import numpy as np
import pandas as pd
import sympy as sy
import matplotlib.pyplot as plt
import timeit

start = timeit.default_timer()

# 读取目标表格csv文件，并用data代表读取到的表格数据
fname = "C:/Users/yc_scorpio/Desktop/test/sea2020.csv"
print("载入数据文件sea2020.csv...")
data = pd.read_csv(fname, encoding='gbk')
m, n = data.shape  # 求数据data的规模 m*n
# 定义探测点距离参数distance，以及海床深度参数deep
distance0 = data.loc[0:(m - 1), ['# 海底海床深度数据']]
distance = np.array(distance0)
distance = [float(x[0]) for x in distance]

deep0 = data.loc[0:(m - 1), ['单位为米']]
deep = np.array(deep0)
deep = [float(x[0]) for x in deep]

# 计算差商y[x0,x1,...,xn]，编程参考教材P91算法DD(x,y,n)
def DD(x,y,n):
    d0 = np.zeros(n+1)
    d0[0] = y[0]
    for k in range(1, n+1):
        for i in range(k, n+1):
            d0[i] = (y[i] - y[i - 1]) / (x[i] - x[i - k])
        for j in range(k, n+1):
            y[j] = d0[j]
    return d0

# d[i]均为二阶差商，有n+1个数据点，len(d)=n+1，返回值为三弯矩方程右端向量d
def cal_d(x,y,n):
    print("二阶差商计算右端向量d[i]...")
    d = np.zeros(n+1)
    # d1=0,dn=0,参考P109公式（3-51）
    d[0] = 0
    d[n] = 0
    m1 = np.zeros(3)
    m2 = np.zeros(3)
    for i in range(1, n):
        m1[0] = x[i-1]
        m1[1] = x[i]
        m1[2] = x[i+1]
        m2[0] = y[i - 1]
        m2[1] = y[i]
        m2[2] = y[i + 1]
        p = DD(m1, m2, 2)
        d[i] = 6*p[2]
    return d

# 定义三弯矩方程系数矩阵aa[i][j]
def cal_aa(n):
    # lam,miu定义均根据P107公式（3-44），lam[0]和miu[n]取值根据公式（3-51）
    print("计算三弯矩方程的系数矩阵a[i][j]...")
    lam = np.zeros(n+1)
    lam[0] = 0
    miu = np.zeros(n+1)
    miu[n] = 0
    # h[i]=100为两标测点距离，n+1为数据点个数
    h = 100*np.ones(n+1)
    for i in range(1, n):
        lam[i] = h[i+1]/(h[i]+h[i+1])
        miu[i] = 1-lam[i]
    aa = np.identity(n+1)*2
    for k in range(n):
        aa[k][k+1] = lam[k]
    for j in range(1, n+1):
        aa[j][j-1] = miu[j]
    return aa

# 高斯消去法解三弯矩方程，求M[i]，编程参考教材P33~34算法GAUSSP
def gaussin1(a,b,n,p,q):
    ## 系数矩阵a，右端向量b，n阶矩阵，上带宽q,下带宽p
    print("解方程得M[j]...")
    cout = 0  # 定义计算次数
    for k in range(n - 1):  # k表示第一层循环，(1，n-1)行
        # 限制条件
        if((k+p)<n-1):
            mm = k+p+1
        else:
            mm = n
        if((k+q)<n):
            nn = k+q+1
        else:
            nn = n
        if (a[k][k] == 0):
            print("no answer")
        for i in range(k + 1, mm):  # i表示第二层循环,(k+1,n)行,计算该行消元的系数
            l = a[i][k] / a[k][k]  # 计算
            cout += 1
            for j in range(k, nn):  # 第i行进行消去
                a[i][j] -= l * a[k][j]
                cout += 1
            b[i]-= l * b[k]
    # 回代求出方程解
    x = np.zeros(n)
    x[n - 1] = b[n - 1] / a[n - 1][n - 1]  # 先算最后一位的x解
    for i in range(n - 2, -1, -1):  # 依次回代倒着算每一个解
        if((i+q)<n):
            nn = i+q+1
        else:
            nn = n
        for j in range(i + 1, nn):
            b[i] -= a[i][j] * x[j]
            cout += 1
        x[i] = b[i] / a[i][i]
        cout += 1
    print("显示前20位数据点二阶导数M[i]...")
    for i in range(20):  # 显示前20位计算解，保留小数点后三位
        print("M[" + str(i + 1) + "] = ", np.around(x[i],3))
    print("解方程浮点计算次数", "≈", cout)
    return x
# 计算每一段(x[i-1],x[i])的三次函数式s_x(n)，为曲线积分做准备
def cal_fi(i):
    x = np.zeros(51)
    y = np.zeros(51)
    for j in range(51):
        x[j] = distance[j]
        y[j] = deep[j]
    h = np.zeros(51)
    h[i] = x[i]-x[i-1]  # 相邻两点间距
    # 每一段三次函数定义为s_x(n)，p为自变量，参考P107公式（3-41）
    s_x = M[i-1]*((x[i]-p)**3)/(6*h[i])+M[i]*((p-x[i-1])**3)/(6*h[i])+\
          (y[i-1]-M[i-1]*(h[i]**2)/6)*(x[i]-p)/h[i]+(y[i]-M[i]*(h[i]**2)/6)*(p-x[i-1])/h[i]
    return s_x

# 定义画图函数draw()，i代表第i个区间x[i]~x[i-1]，并求出曲线长度l1
def draw():
    ## 开始画图
    print("画图中...")
    area_list = []  # 存储每一微小步长的曲线长度
    l1 = 0
    for i in range(1, len(distance)):
        p1 = cal_fi(i)  # 求出第i段三次函数
        # 根据公式p1画出每一段曲线
        # 并求出每一段的曲线弧长l[i]，l[i]存储于列表area_list
        xx = np.arange(0, 100, 0.1)+100*(i-1)
        yy = np.zeros(len(xx))
        yy[0] = p1.evalf(subs={p: xx[0]})
        for j in range(1, len(xx)):
            yy[j] = p1.evalf(subs={p: xx[j]})
            l = (yy[j] - yy[j - 1]) ** 2 + 0.1 ** 2
            l = np.sqrt(l)
            area_list.append(l)
        plt.plot(xx, yy)
        l1 = l1 + sum(area_list)
        area_list = []
    print("画图结束...")
    print("计算曲线长度 = ", np.around(l1, 1), "m")

# 画布参数设置
def draw_setting():
    plt.scatter(distance, deep, s=10)  # 画出原始数据散点图
    # 默认英文，如果添加了中文标题，则会输出一堆乱码，解决方法就是加上下面这两行代码
    plt.rcParams['font.sans-serif'] = ['SimHei']
    plt.rcParams['axes.unicode_minus'] = False
    plt.xlabel("distance/m", fontsize=14)
    plt.ylabel("deep/m", fontsize=14)
    plt.title("三次样条插值曲线", fontsize=16)
    # plt.legend(['curve-fitting'])
    plt.ylim(350, 300)
    plt.xlim(-100, 5100)
    plt.show()

if __name__ == '__main__':
    d = cal_d(distance, deep, len(distance)-1)  # 三弯矩方程右端向量d[]
    aa = cal_aa(len(distance)-1)  # 三弯矩方程系数矩阵aa[][]
    M = gaussin1(aa, d, 51, 1, 1)  # 高斯消去法解方程得M[]，并输出前20个数据点的二阶导数
    p = sy.Symbol('p')  # 将p定义为自变量
    draw()   # 画图，并求出曲线长度
    draw_setting()  # 画图参数设置
    end = timeit.default_timer()
    tmp = np.around(1000 * float(end - start), 2)
    print("运行时长 =", tmp, "ms")
